# Chargement des packages
library(dplyr)
library(sf)
library(spdep)
library(RColorBrewer)
# Import des données
iris<-st_read("./fonds/iris_franceentiere_2021/iris_franceentiere_2021.shp",
options="ENCODING=WINDOWS-1252")
## options: ENCODING=WINDOWS-1252
## Reading layer `iris_franceentiere_2021' from data source
## `/home/onyxia/work/Stat_Spatiale_TP6/TP6_Stat_spatiale/fonds/iris_franceentiere_2021/iris_franceentiere_2021.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 16351 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -63.15332 ymin: -21.38963 xmax: 55.83665 ymax: 51.0668
## Geodetic CRS: WGS 84
data<-read.csv2("./data/BASE_TD_FILO_DISP_IRIS_2018.csv",sep=";",dec=".")
# Jointure
marseille<-iris %>%
filter(substr(depcom,1,3)=="132") %>%
left_join(
data %>% select(code=IRIS,DISP_MED18),
by="code"
)
## Classes 'sf' and 'data.frame': 393 obs. of 5 variables:
## $ code : chr "132010101" "132010102" "132010103" "132010104" ...
## $ libelle : chr "La Bourse" "Thubaneau" "Colbert-Providence" "Bernard du Bois" ...
## $ depcom : chr "13201" "13201" "13201" "13201" ...
## $ DISP_MED18: int 16270 11320 11470 11920 14110 20030 15990 NA 11750 13140 ...
## $ geometry :sfc_MULTIPOLYGON of length 393; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:57, 1:2] 5.37 5.37 5.37 5.37 5.37 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
## ..- attr(*, "names")= chr [1:4] "code" "libelle" "depcom" "DISP_MED18"
marseille <- marseille %>%
st_transform(2154)
summary(marseille$DISP_MED18)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 10000 16115 19255 19767 23490 38500 57
boxplot(marseille$DISP_MED18)
hist(marseille$DISP_MED18, breaks = seq(10,40,1)*1e3) #par iris
# par arrondissement
boxplot(marseille$DISP_MED18 ~marseille$depcom)
#On peut faire un test de Fisher pour analyser la variance
#variance intra/interclasse. Est-ce que la Variance inter est plus forte que la variance intra ?
summary(aov(DISP_MED18 ~ depcom, data = marseille))
## Df Sum Sq Mean Sq F value Pr(>F)
## depcom 15 5.759e+09 383941516 30.7 <2e-16 ***
## Residuals 320 4.003e+09 12508008
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 57 observations deleted due to missingness
# <2e-16
#on rejette l'hypothèse nulle car la p-valeur est très faible. Il y a une différence entre les arrondissements
library(ggplot2)
ggplot(marseille) +
geom_boxplot(aes(x=DISP_MED18, y = depcom)) +
geom_vline(xintercept = mean(marseille$DISP_MED18), linetype = "dashed", size = 1, col = "darkred", alpha= 0.5)+
labs(y="arrondissements")+
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## Warning: Removed 57 rows containing non-finite values (`stat_boxplot()`).
## Warning: Removed 1 rows containing missing values (`geom_vline()`).
marseille <- marseille %>%
filter(!is.na(DISP_MED18))
# tidyr::drop_na() # autre façon de retirer les valeurs manquantes
plot(marseille["DISP_MED18"]) #analyse continue
mapview(
marseille,
z= c("DISP_MED18"),
layer.name = "Revenu_median",
alpha.regions = 0.35, #transparence de la couche
label = "code"
)
plot(marseille["DISP_MED18"], breaks = "quantile", nbreaks = 10) #analyse avec discrétisation automatique
plot(marseille["DISP_MED18"], breaks = "jenks") #analyse avec discrétisation automatique
sample() et à partir de la variable
DISP_MED18 du fond des iris de Marseille. Vous stockerez ce
vecteur dans une nouvelle variable du fond des iris nommée
DISP_MED18_ALEA.set.seed(1793)
marseille <- marseille %>% mutate(DISP_MED18_ALEA = sample(DISP_MED18))
library(gridExtra)
##
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
##
## combine
map1 <- ggplot(marseille) +
geom_sf(aes(fill = DISP_MED18, col = DISP_MED18)) +
guides(col = "none") +
scale_fill_viridis_c("ORIG") +
scale_color_viridis_c("") +
theme_void()
map2 <- ggplot(marseille) +
geom_sf(aes(fill = DISP_MED18_ALEA, col = DISP_MED18_ALEA)) +
guides(col = "none") +
scale_fill_viridis_c("ALEA") +
scale_color_viridis_c("") +
theme_void()
# on met en forme les deux cartes
grid.arrange(map1,map2, layout_matrix = matrix(1:2, nrow=2))
la comparaison des deux cartes semble suggérer que la carte représentant la distribution réelle des revenus est très différente de la carte d’une distribution aléatoire. LE phénomène semble spatialement corrélé. 6. Pour corroborer la conclusion du 5., nous allons mesurer et tester le phénomène d’autocorrélation spatiale.
Un phénomène est autocorrélé spatialement quand la valeur de la variable étudiée à un endroit donné est plus liée aux valeurs de ses voisins plutôt qu’à celles des autres. On parle d’autocorrélation positive si des voisins ont tendance à prendre des valeurs similaires et d’autocorrélation négative si des voisins ont tendance à prendre des valeurs différentes.
Quel type d’autocorrélation spatiale, le phénomène étudié semble-t-il avoir ?
Pour étudier le phénomène, il nous faut construire une matrice de voisinage. Il existe plusieurs façons de définir le voisinage de nos iris. Dans un premier temps, nous allons définir le voisinage par la contiguïté : deux iris sont voisins s’ils sont contigus.
Pour limiter la taille des objets créés, nous allons travailler avec des listes plutôt qu’avec des matrices carrées.
Extraire la liste des voisins de chaque Iris. Pour cela, vous
utiliserez la fonction spdep::poly2nb(). Par défaut, il
s’agit de la contiguité dite QUEEN qui reprend les
mouvements de la Reine aux échecs. Prenez connaissance de l’objet créé
et réalisez un résumé de l’objet en sortie avec la fonction
summary().
Créez une liste de poids à partir de la liste de voisins
précédemment créée. Pour cela, utilisez la fonction
spdep::nb2listw(), avec l’argument
zero.policy=TRUE pour intégrer les Iris n’ayant
potentiellement pas de voisins (par défaut, la fonction exclut les
observations sans voisin).
Prenez connaissance de l’objet créé avec la fonction
str() et l’argument max.level = 1 et réalisant
un résumé de la liste avec la fonction summary().
Vérifiez que la somme des pondérations associées à chaque pondération est égale à 1.
Créer une variable des revenus disponibles centrés réduits avec
la fonction scale(). Vous la nommerez
DISP_MED18_STD et l’ajouterez au fond des iris de
Marseille. Vous vérifierez que la nouvelle variable est bien centrée
(moyenne = 0) et réduite (SD = 1).
Dresser le diagramme de Moran avec la fonction
moran.plot() à partir de la variable standardisée (utiliser
la fonction as.numeric() si un problème apparaît). Le
second argument à préciser (listw) correspond à la liste
des poids des voisins que vous avez créée précédemment.
Le diagramme de Moran représente, pour chaque observation (ici un
Iris), croise deux informations :
Interprétez les quatre cadrans du diagramme.
Calculez cet indice et sa significativité avec la fonction
spdep::moran.test() utilisée de la façon suivante :
moran.test(marseille$DISP_MED18_STD, ponderation, randomisation = TRUE).
Le dernier argument signifie que la distribution observée est comparée à
une distribution aléatoire obtenue par permutation des valeurs
observées.
Interpértez le résultat obtenu : confirme-t-il ou non votre hypothèse ?
L’indice de Moran est un indicateur global de mesure de l’autocorrélation. Mais, ce phénomène peut connaître une intensité très différente localement. Dans certains endroits de la ville de Marseille, la ressemblance des voisins peut être très forte et à d’autres endroits plus lâche. Des indicateurs locaux de mesure de l’autocorrélation spatiale sont nécessaires pour compléter l’analyse de la distribution spatiale des revenus disponibles médians à Marseille. Nous calculerons pour cela les LISA (Local Indicators of Spatial Association), ou I de Moran locaux.
Calculez les Lisa avec la fonction
spdep::localmoran() et stockez le résultat dans un objet
appelé mars_rev_lisa.
Etudiez l’objet obtenu, en utilisant notamment les fonctions
class(), str(.,max.level=1) et
summary().
Quelle est la moyenne des indicateurs locaux (\(I_i\))?
L’interprétation d’un Lisa est tout à fait similaire à l’indice global. Si l’indicateur local d’un Iris donné est positif, cela signifie qu’il est entouré d’Iris ayant des niveaux de revenus similaires. S’il est négatif, cela indique qu’il est plutôt entouré d’Iris ayant des niveaux de revenus différents (opposés).
Combien d’indicateurs locaux sont-ils négatifs ?
Nous cherchons à représenter les Lisa sur une carte des Iris
marseillais. Pour cela, ajouter les Lisa comme une nouvelle variable du
fond des iris, variable que vous nommerez LISA.
Interprétez ce que vous voyez sur la carte.
Comme pour le I de Moran, il est nécessaire avant d’aller plus
loin dans l’interprétation, de savoir si les Lisa calculés sont
significativement différents de zéro. Dans l’objet
mars_rev_lisa, repérez la colonne correspondant à la
pvaleur du test associé à la mesure des Lisa et placez-la dans une
nouvelle variable du fond des Iris intitulée
LISA_PVAL.
Combien de LISA sont-ils significativement différents de zéro pour un niveau de confiance à 95% ?
Représentez sur une carte la p-valeur des LISA en choisissant les bornes d’intervalles suivantes : 0,0.01,0.05,0.1,1.
Les zones précédemment repérées sur la carte des LISA font-elles parties des zones où les LISA sont les plus significatifs ?